
% This code produces two plots showing the catchment areas where 9 plants
% are located in a square for two cases of constant b_0 and varied b_o.
% Last update 03/01/2023
clear all
close all
clc

figdir = '..\..\Plots\';

points = 500 ;

g_x = linspace(0.0,1.0,points) ;
g_y = linspace(0.0,1.0,points) ;

d = @(x) exp(((x(1,1)-x(3,1)).^2 + (x(2,1)-x(4,1)).^2).^(1/2)) ;

p1 = [1/6 5/6] ;
p2 = [1/2 5/6] ;
p3 = [5/6 5/6] ;
p4 = [1/6 1/2] ;
p5 = [1/2 1/2] ;
p6 = [5/6 1/2] ;
p7 = [1/6 1/6] ;
p8 = [1/2 1/6] ;
p9 = [5/6 1/6] ;

%% Figure 1
weights = [1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0] ;

for i=1:points
    for j=1:points
        x1 = [p1(1) p1(2) g_x(i) g_y(j)]' ;
        x2 = [p2(1) p2(2) g_x(i) g_y(j)]' ;
        x3 = [p3(1) p3(2) g_x(i) g_y(j)]' ;
        x4 = [p4(1) p4(2) g_x(i) g_y(j)]' ;
        x5 = [p5(1) p5(2) g_x(i) g_y(j)]' ;
        x6 = [p6(1) p6(2) g_x(i) g_y(j)]' ;
        x7 = [p7(1) p7(2) g_x(i) g_y(j)]' ;
        x8 = [p8(1) p8(2) g_x(i) g_y(j)]' ;
        x9 = [p9(1) p9(2) g_x(i) g_y(j)]' ;
        dist = [d(x1) d(x2) d(x3) d(x4) d(x5) d(x6) d(x7) d(x8) d(x9)]./weights ;
        [min_dist c(i,j)] = min(dist) ;
        clear dist min_dist
    end
end

x_loc = [x1(1) x2(1) x3(1) x4(1) x5(1) x6(1) x7(1) x8(1) x9(1)] ;
y_loc = [x1(2) x2(2) x3(2) x4(2) x5(2) x6(2) x7(2) x8(2) x9(2)] ;

for x=2:(points-1)
    for y=2:(points-1)
        if (c(x,y)~=c(x,y+1)) & (c(x,y+1)~=0)
            c(x,y)=0;
        elseif (c(x,y)~=c(x,y-1)) & (c(x,y-1)~=0)
            c(x,y)=0;
        elseif (c(x,y)~=c(x+1,y)) & (c(x+1,y)~=0)
            c(x,y)=0;
        elseif (c(x,y)~=c(x-1,y)) & (c(x-1,y)~=0)
            c(x,y)=0;
        end
    end
end

for x=1:points
    for y=1:points
        if c(x,y)~=0
            c(x,y)=1;
        end
    end
end
fig1 = figure(1);
h=pcolor(g_x,g_y,c);
set(h,'LineStyle','none')
hold on
scatter(x_loc,y_loc,'filled','k')
mymap = [0 0 0 ; 1 1 1];
colormap(mymap)
labels = string(weights)
for num = 1:length(labels)
    labels(num)=strcat('(', labels(num), ')')
end
set(gca,'xtick',[])
set(gca,'xticklabel',[])
set(gca,'ytick',[])
set(gca,'yticklabel',[])
text(y_loc,x_loc,labels,'VerticalAlignment','top','HorizontalAlignment','center')
rectangle('Position',[0 0 1 1])


%saveas(fig1,[figdir 'Figure4_left'],'png');


%% Figure 2


points = 500 ;

g_x = linspace(0.0,1.0,points) ;
g_y = linspace(0.0,1.0,points) ;

d = @(x) exp(((x(1,1)-x(3,1)).^2 + (x(2,1)-x(4,1)).^2).^(1/2)) ;

p1 = [1/6 5/6] ;
p2 = [1/2 5/6] ;
p3 = [5/6 5/6] ;
p4 = [1/6 1/2] ;
p5 = [1/2 1/2] ;
p6 = [5/6 1/2] ;
p7 = [1/6 1/6] ;
p8 = [1/2 1/6] ;
p9 = [5/6 1/6] ;

weights = [1.2 1.0 1.0 1.0 1.2 1.0 1.0 1.0 0.85] ;


for i=1:points
    for j=1:points
        x1 = [p1(1) p1(2) g_x(i) g_y(j)]' ;
        x2 = [p2(1) p2(2) g_x(i) g_y(j)]' ;
        x3 = [p3(1) p3(2) g_x(i) g_y(j)]' ;
        x4 = [p4(1) p4(2) g_x(i) g_y(j)]' ;
        x5 = [p5(1) p5(2) g_x(i) g_y(j)]' ;
        x6 = [p6(1) p6(2) g_x(i) g_y(j)]' ;
        x7 = [p7(1) p7(2) g_x(i) g_y(j)]' ;
        x8 = [p8(1) p8(2) g_x(i) g_y(j)]' ;
        x9 = [p9(1) p9(2) g_x(i) g_y(j)]' ;
        dist = [d(x1) d(x2) d(x3) d(x4) d(x5) d(x6) d(x7) d(x8) d(x9)]./weights ;
        [min_dist c(i,j)] = min(dist) ;
        clear dist min_dist
    end
end

x_loc = [x1(1) x2(1) x3(1) x4(1) x5(1) x6(1) x7(1) x8(1) x9(1)] ;
y_loc = [x1(2) x2(2) x3(2) x4(2) x5(2) x6(2) x7(2) x8(2) x9(2)] ;

for x=2:(points-1)
    for y=2:(points-1)
        if (c(x,y)~=c(x,y+1)) & (c(x,y+1)~=0)
            c(x,y)=0;
        elseif (c(x,y)~=c(x,y-1)) & (c(x,y-1)~=0)
            c(x,y)=0;
        elseif (c(x,y)~=c(x+1,y)) & (c(x+1,y)~=0)
            c(x,y)=0;
        elseif (c(x,y)~=c(x-1,y)) & (c(x-1,y)~=0)
            c(x,y)=0;
        end
    end
end

for x=1:points
    for y=1:points
        if c(x,y)~=0
            c(x,y)=1;
        end
    end
end

fig2 = figure(2);
h=pcolor(g_x,g_y,c);
set(h,'LineStyle','none')
hold on
scatter(x_loc,y_loc,'filled','k')
mymap = [0 0 0 ; 1 1 1];
colormap(mymap)
labels = string(weights)
for num = 1:length(labels)
    labels(num)=strcat('(', labels(num), ')')
end
set(gca,'xtick',[])
set(gca,'xticklabel',[])
set(gca,'ytick',[])
set(gca,'yticklabel',[])
text(y_loc,x_loc,labels,'VerticalAlignment','top','HorizontalAlignment','center')
rectangle('Position',[0 0 1 1])

%saveas(fig2,[figdir 'Figure4_right'],'png');
